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We compute the electrostatic potential at the surface, or zeta potential C, of a 
charged particle embedded in a colloidal suspension using a hybrid mesoscopic 
model. We show that for weakly perturbing electric fields, the value of ^ obtained 
at steady state during electrophoresis is statistically indistinguishable from C, in 
thermodynamic equilibrium. We quantify the effect of counterions concentration 
on C,. We also evaluate the relevance of the lattice resolution for the calculation of 
C, and discuss how to identify the effective electrostatic radius. 
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1. Introduction 

A quantitative understanding of the electrostatic interactions between charged 
macroions in solution is fundamental to comprehend a plethora of physical phe- 
nomena spanning from biology to material science, for example the macroscopic and 
rheological properties of colloidal suspensions[l,2]. The details of such interactions 
are difficult to capture, as the effective interactions are determined by the interplay 
between the different components dissolved in solution (macroions, counterions, salt 
ions) and the solvent dielectric response[3]. In addition, when the system is driven 
out of equilibrium, for example by an external electric field causing electrophoretic 
flow, hydrodynamic interactions between the solvent and solute species must also 
be included and can alter the equilibrium electrostatic interactions. 

The simplest yet very useful model to describe electrolyte solutions is the Poisson- 
Boltzmann (PB) approach[l,4], which builds on a continuous description of the 
electrolyte, characterized in terms of the anion and cation local densities. Within 
this framework, it is possible to add charged macroscopic objects, with fixed surface 
charges, by accordingly changing the electrostatic boundary conditions of the PB 
equations. Notwithstanding PB is a mean-field theory that defines ions as point-like 
and therefore neglects excluded volume effects and spatial correlation due to their 
finite size, it has been successfully adopted to describe various physical systems, for 
example to derive the electrostatic part of DLVO interparticle potential[l] which 
explains the interactions of weakly charged particles in solution at low volume frac- 
tion. 

The zeta potential, C, defined as the electrostatic potential at the colloid sur- 
face, or slippling plane, plays a central role for charged colloidal dispersions, as it 
indirectly provides an estimate of the magnitude of the electric field between parti- 
cles that results from the combined effect of ionized charged groups sitting on the 
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particle, counterions released into solution and salt ions. Theoretically, Olishima 
et a/. [5] obtained an exact analytic expression relating the surface charge density 
to the zeta potential for an infinitely dilute spherical colloidal suspension by care- 
fully approximating the non-linear PB equation. This expression is of limited use, 
since the infinite dilution regime is difficult to achieve. Moreover, experimentally 
C, is normally derived out of equilibrium, from electrophoretic mobility measure- 
ments. Since, a priori, the values of C, in equilibrium and at steady state will differ, 
an accurate description of hydrodynamics must be included in order to derive the 
relationship between C and the colloid mobility. 

Analytical results for the electrophoretic flow at finite dilution cannot in gen- 
eral be obtained. Computationally, different models have been used to describe 
electrophoresis[6-9,ll]. Lobaskin et aZ.[8] and Diinweg et a/[9], using a lattice Boltz- 
mann[10] (LB) solver coupled to molecular dynamics of ions, calculate the particle 
mobility for low to zero salt concentration explicitly accounting for finite ions size, 
but without providing results for C,. Kim et al\l\\ show agreement with Ohshima's 
results, but only for weakly charged particles at low salt concentration. Cell mod- 
els[12] are also employed to describe electrokinetics and to derive mobility versus 
C, curves; however, they rely on somewhat arbitrary boundary conditions both for 
electrostatics and hydrodynamics. 

In this paper, we compute C, using a mesoscopic model that couples a Navier- 
Stokes solver to the convection and diffusion of ions on a lattice. We will first 
investigate the difference between C, values at thermodynamic equilibrium and at 
electrophoretic steady state. We then accurately analyze the dependence of C, on 
colloid volume fraction, $, and salt and counterions concentration, and compare 
them with predictions from PB at infinite dilution derived by Ohshima et al.\b\. 
Finally, we will address how to identify the electrostatic radius of a colloid on 
a lattice, as it has been done previously for the hydrodynamics radius[13]. In the 
following section, we briefly introduce our model, the reference equations and results 
that can be obtained using PB equation for infinite dilution and at finite volume 
fractions. We present in section [3] results for C as a function of the volume fraction, 
salt concentration and lattice refinement. Finally, we draw our main conclusions in 
the last section. 



In order to derive a correct relation between the mobility and the zeta potential 
at steady state during electrophoresis, our model includes hydrodynamic interac- 
tions and electric forces due to local charge densities. 



2. Model 



(a) Electrophoretic flow 



dpk 
dt 



(2.1) 




(2.2) 



k 




(2.3) 



k 
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where Dk-,Zk,k = +, — are the diffusivities and valences of positive and negative 
ions, p,v,pid and r/ correspond to the solvent density, velocity, ideal pressure and 
shear viscosity, e is the electron charge, jS'"^ = fc^T is the Boltzmann factor, ip 
is the electrostatic potential and e the solvent permittivity. Eq. 12.11 expresses ions 
mass conservation during diffusion and advection, coupling ion dynamics to solvent 
motion. This is in turn described by the Navier-Stokes equation for a viscous fluid 
which couples solvent dynamics to electrostatic forces due to local charge density, 
Eq. 12.21 Finally, the Poisson equation enforces the electrostatic coupling between 
the charged species and the embedded solid particles. We solve these electrokinetics 
equations combining an LB approach for the momentum dynamics with a numeri- 
cal solver for the discretized convection-diffusion dynamics, Eq. (j2.ip . based on link 
fluxes[14]. Each lattice site is labeled as solid or fluid and we consider here a single 
spherical object of radius r = a embedded in a cubic lattice of volume with pe- 
riodic boundary conditions. Suspended particles are therefore resolved and interact 
with the neighbouring fluid through bounce-back[13]. The electrostatic potential is 
computed using a successive over- relaxation scheme (S0R)[16]. This method has 
been successfully applied to analyze different dynamical processes involving suspen- 
sions of charged objects[15,17,18,19]. 



(6) Poisson- Boltzmann electrostatics 

At equilibrium, using a mean-field approach that neglects excluded volume ef- 
fects and correlation between ions, PB equation reads[4] 

„9 SirezpQ . , ,ezip(r). 

VV r- = ^ sinh , 2.4 

where a symmetric electrolyte z+ = z_ = z has been used for simplicity and po, 
the uniform macroscopic counterion and colon number concentration far from the 
colloid p+ = P- = POj is assumed to be equal to that of an electrolyte reservoir 
with dissolved salt Cgait = 2/9o- 

An analytical solution of eq. (|2.4p for a single particle is not available, as the 
PB equation can be solved analytically only for a few symmetrical configurations. 
Ohshima et al.[5] analyzed eq. (j2.4p around a spherical colloid of charge Q using a 
perturbative approach and derived a quasi-exact analytical expression for C, 

nog(cosh2(c74))' 



Q = 2t-^ sinh(C/2) 

Acts 



1 



(Ka)cosh^(C/4) (Ka)2sinh^(C/2) 



(2.5) 



where n — \ = '^n^ i'' inverse of the Debye length, l^ — . \ rp the 



Bjerrun length and ^ = ^^^^ adimensional zeta potential. 

Useful global information can be obtained using the Debye-Hiickel approxima- 
tion, i.e. linearizing eq. ^T^. When < 1.0, 

V^^(r) = lq!^^(f) = (2.6) 

from which one can derive the electrostatic potential around a charged spherical 
colloid 

<p(r) = ^cxp(-('"-"'/^°) . (2.7) 
r 
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In the linear regime, a charged particle in the presence of salt develops a screened, 
Yukawa-like electrostatic potential, with a characteristic length scale oi \d and 
strength given by the zeta potential, C. 

Ea. (j2.4p is strictly valid when the amounts of positive and negative charges 
dissolved in solution are equal, i.e. when the counterions concentration is negligible. 
Experimentally, this can be obtained separating the colloidal particles from a bulk 
electrolyte solution with a semi-permeable membrane. For the experimental set-ups 
with non-zero concentration of counterions, eq. (j2.6p becomes 



7y — + (1 + 1^ — mn 



(2.8) 



where pso, Pco are the reference salt and counterion concentrations in the reservoir. 
Eq. (|2.8p stresses the fact that, especially at low salt concentration Cs (low pso) 
and high volume fraction (j) ~ (high Pco), a deviation from the theoretical 

Debye-Hiickel results can be observed. 

Our aim is to measure ^ for different experimental regimes, analyze its behavior 
beyond the linear approximation (when -j^^ > 1) and understand how to overcome 
the intrinsic inaccuracies in the location of the electrostatic radius of any lattice 
model. In the LB approach, we will refer to ( as the average electrostatic potential 
between lattice points pairs (boundary links) I that link a colloid and a fluid node 
calculated half-way between the two nodes using linear approximation. 

and will analyze when such an assumption is representative of the colloid 



3. Results 

We have performed simulations for a single spherical particle of radius a embedded 
in cubic boxes with periodic boundary conditions, also varying the concentration 
of added salt. In order to avoid finite-size effects, for a given salt concentration we 
use a lattice length L at least three times the corresponding Debye length, Xu. We 
have distributed uniformly the colloid charge Q on the lattice sites the colloid fills, 
and choose Q to correspond to the predicted value of Oshima, according to eq. (j2.5p . 
We will analyze both a magnitude of Q which lies in the linear regime, Co = 1-0, 
and one for which the colloid is within the nonlinear regime, Co = 5.0. We finally 
set the external electric field to E ~ 0.01, well within the linear regime E 

The values for ( obtained from electrophoresis are in principle different from 
those computed in equilibrium. This is analogous to the different measurements 
for the size of a solvated particle, at equilibrium or at hydrodynamic steady state 
(hydrodynamic radius). We have therefore computed Cmp at thermodynamic equi- 
librium, i.e. with equilibrated charges and fluid at rest and dynamically, when a 
steady state due to the external electric field is reached. In all simulations, our 
results show that the differences between the two values of ( are smaller than the 
statistical precision. This similarity is expected because of the small magnitude of 
the perturbing electric field and colloid Peclet number, Pe ~ , where Vcou is 
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Figure 1. Left(Right): Electrostatic potential (mp calculated halfway between boundary 
links (see text for definition) using a low(higli)-resolution lattice. Colloidal charges for 
different salt concentrations {ka = 0.5, 2.25, 4.5) have been assigned so that the resulting 
theoretical zeta potential (eq. H2.5[l ) is within the linear electrostatic, Debye-Hiickel regime 
(Co ~ 1.0, dashed line) 



the velocity of the colloid at steady state. For Pe > 1 the deformation of the salt 
cloud around the colloid induces significant deviations in the electrophoretic mo- 
bility [18]; we can hence envisage a corresponding departure of the dynamic C from 
its equilibrium counterpart. Typically, as for simulations here, Pe <^ 1 and particle 
mobility does not depend on Pe, hence the measured ^ values at equilibrium and 
steady state are statistically indistinguishable. 

As the only theoretical result available, Co, is valid for a single particle in equi- 
librium at infinite dilution, we will compare the simulation results with eq. (j2.5p to 
assess the role of the counterion concentration at finite volume fraction. In addi- 
tion, in order to analyze to which extent the midpoint stands as a proper definition 
for the electrostatic radius, we run a pair of simulations for each colloid charge, 
volume fraction and salt concentration used, changing the resolution of the colloid 
in the lattice from low, using colloid radius a = 4.5, to high a = 9.0, also doubling 
the corresponding box sizes (e.g. from L = 30 at low resolution to L = 60 at high 
resolution for the highest volume fraction). 

FigH] displays (mp for a weakly charged colloid (Co = 1-0) as a function of 
the system size and salt concentrations, where the left(right) panel corresponds to 
a low(high) colloid resolution. The top panels show that the deviation of C from 
Oshima's prediction varying $ appears indirectly only through the corresponding 
change in the co- and counter-ion concentrations because the deviation decreases 
when the salt concentration is increased. Therefore, according to eq. (j2.8p . we can 
gain understanding by analyzing the dependence of Cmp on the ratio between the 
concentration of counterions released by the particle, Cdons (which decreases as $ 
increases ) and the concentration of salt Csait dissolved in the system (which is 
independent of $). The bottom panels confirm an agreement within 10% between 
Cmp and Co for ^ 0.1. In addition, comparing the left and right panels, we 

conclude that in the linear regime the lattice resolution does not play a significant 
role, as for a given salt concentration Cmp does not differ significantly and the use 
of a more computationally expensive refined lattice only reduces the error bars. In 
addition, we note from bottom-rigth panel that discretization effects become more 
relevant when increasing the salt concentration. This sensitivity is associated to 
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Figure 2. Left(Right): Electrostatic potential C,mp calculated halfway between boundary 
links (see text for definition) using a low(high)-resolution lattice. Colloidal charges for 
different salt concentrations {ka = 0.5, 2.25, 4.5) have been assigned to obtain a theoretical 
zeta potential (eq. (|2.5p . dashed line) C,o = 5.0, within the nonlinear regime 



the reduction of Ad and the corresponding loss of resolution in the electrostatic 
potential around the colloid. 

Fig. [5J organized analogously to Fig. [1] shows results for a strongly charged 
colloid, C^o = 5.0, for which the linearized Debye-Hiickel approximation does not 
hold. We observe that the dependence of Cmp on $ enters again indirectly through 
the relative changes in Cdonsi and that the convergence toward Co can only be 
expected for high Csait- However, as opposed to the linearized regime, adding salt 
can lead up to a 20% overestimation of Cmp for the highest salt concentration 
analyzed {ka = 4.5), especially when using a low resolution lattice (left panels). 
In this strong coupling regime the electrostatic potential decays faster than A£)[l]. 
This strong nonlinear behavior invalidates the midpoint as the natural choice for 
the electrostatic radius, which develops a dependence on the system parameters 
even when the effects of finite $ and Cdons are negligible. 

Fig. [3] displays the measured electrostatic potential at the midpoint between a 
solid and fluid node, Ca/p, together with the two extreme situations where the elec- 
trostatic potential is averaged over the solid and the fluid nodes corresponding to 
the boundary links. One can observe, as expected, that Cmp always lies in between 
the other two estimates of C. The differences observed between the left and right 
panels indicate that the increase in resolution does not affect significantly Ca/p, 
while decreases the inaccuracy in the estimate of the limiting values for C,. Gener- 
ically, the better resolved the colloidal particle (and the corresponding decay of 
the electrostatic potential, i^) the less spreading between the different electrostatic 
potential estimates. When we increase Cdons, Oshima's prediction does no longer 
hold and eq. (|2.5p cannot be used to identify the electrostatic radius as can be 
appreciated in the rightmost values of C in the bottom panels. However, the weak 
dependence of tp in the boundary fluid nodes suggests that when Cmp is no longer 
valid, we can still identify a electrostatic radius slightly larger than the one associ- 
ated to Cmp- Analogously to the hydrodynamic radius, an ad hoc calibration of the 
electrostatic radius at all salt concentrations must be done to perform quantitative 
studies of the electrokinetics of colloidal suspensions with LB. 
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Figure 3. Left[Right]: Average electrostatic potential calculated using low[high] resolution 
lattice at midpoint (squares), solid (triangle-up) and fluid (triangle-down) nodes pertaining 
to boundary links in all simulation set-ups (ka = 0.5, 2.25, 4.5, (o = 1-0, 5.0). Top[Bottom] 
panels correspond to the case of a weakly [strongly] charged colloid. 



4. Conclusions 

In this paper, we have analyzed how to determine the electrostatic potential at the 
surface of a particle, or zeta potential C, for a colloidal suspension. To this end, 
we have made use of a hybrid mesoscopic model which couples a discrete lattice 
formulation of Boltzmann's kinetic equation for the solvent to a discrete solution 
of the convection-diffusion equation for the charged ion species dissolved in the 
fluid which implies treating counter- and salt ions as scalar fields at the Poisson- 
Boltzmann (PB) level. We have found that for weakly perturbing electric fields, it 
is not possible to distinguish between the equilibrium and dynamic zeta potentials. 
We expect differences will arise when the deformation of the charge layer around 
the colloid becomes significant, a scenario which can be achieved, e.g. when the 
Peclet number of the dissolved ions is not negligible. 

In particular, the comparison with Oshima's expression for ^ at infinite dilution 
has allowed us to carry out quantitative checks to validate the code performance. 
We have seen that the counterion concentration has a significant effect on (, leading 
to a decrease in its magnitude both in the linear and nonlinear regimes away from 
Oshima's result. We have shown that the ratio between counterion and salt con- 
centration controls the departure from Oshima's prediction, and that its prediction 
works reasonably well for Cdons / Csait ^ 0.1. We have also assessed the relevance of 
the lattice resolution and have quantified its effects. We have seen that the mean of 
( over the boundary nodes which determine the colloidal shape is in general a good 
estimate of the electrostatic radius. However, for highly charged colloids a more 
refined choice of the particle radius will be in general needed. The effective radius 
will be slightly larger than predicted from Cmp due to the nonlinear decay of the 
electrostatic potential around the particle. This effective electrostatic radius, which 
needs to be calibrated as a function of the salt concentration and particle radius, 
will in general differ from the particle hydrodynamic radius and requires a separate 
analysis. The overall dependence of the electrostatic radius both on the applied 
field and ion concentrations is weaker than the one observed for the hydrodynamic 
radius. In most situations we expect that an equilibrium calibration correcting from 
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the finite counterion concentration will provide a quantitative estimate for the elec- 
trokinetics of colloidal suspensions. 
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